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Abstract. 

We describe our methods of fitting a smoothing function to observational range differences 
from a predicted orbit, by deducing corrections to the orbit in the radial and along-track 
directions. The method has been used on observations of a variety of satellites, and using 
predicted orbits computed both by numerical integration using IRVs as starting values and 
analytically from orbital elements. The along-track corrections to the predicted orbit have 
been succesfully used in the form of time biases to improve subsequent predictions, and a sta- 
tistical test has been devised to ensure that the range residuals may be used to form unbiased 
quicklook normal points. 


1. Introduction. 

It was recommended at the Fifth International Laser Ranging Instrumentation Workshop 
held at Herstmonceux, UK, that laser range normal points should be formed by the stations 
shortly after each satellite pass and transmitted as quicklook data (Gaignebet, 1985). The 
final version of the quicklook normal point format was agreed by the SLR Subcommission of 
the CSTG and published in the April 1990 Satellite Laser Ranging Newsletter (Schutz, 1990). 
The process of forming normal points consists of 2 stages; (a) fitting a smoothing function 
to the observational range differences from a predicted orbit and subtracting the function in 
order to form a flat track of residuals, and (b) rejecting noise and outliers and forming the 
mean values of the accepted residuals in bins spread through the pass. The values of the 
smoothing function and the predicted orbit at the epoch of an actual observation nearest 
the mean epoch of the bin are then added to the mean value of the residuals in the bin. The 
resulting normal points axe then virtually independent of both the smoothing function and 
of the predicted orbit. 

In this report we detail the processes adopted to carry out stage (a) of the process, methods 
for rejecting noise and forming unbiased normal point bin means are developed by Sinclair 
(these proceedings) 

2. Development of the Smoothing Function. 

During the observation of a satellite pass, differences from the predicted range are computed 
for all events detected within the range gate and displayed in realtime, true satellite returns 
appearing as a track of correlated values among the randomly distributed noise events. An 
example of such a plot is given in Figure 1, which shows a pass of ERS-1. The rejection of 
gross noise events is carried out visually using the plot, and the subsequent set of satellite 
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range measurements and some noise events axe passed to the next stage of processing. 

At an early stage of the development of the preprocessing software, the decision was taken 
to use the pass-by-pass observations to sdlve for corrections to the predicted orbit, and 
in particular to monitor the along-track error, in order to improve subsequent prediction 
accuracy. Thus the possible use of an arbitrary polynomial for the smoothing function was 
not considered. However, corrections to the predicted orbital elements of the satellite would 
be difficult to apply to orbits computed by numerical integration and using IRVs as starting 
values, and similarly corrections to the IRVs would not be applicable to orbits computed 
analytically. Thus it was decided to express orbital corrections as along-track, across-track 
and radial displacements to the predicted position. Such corrections are then readily applied 
to the geocentric orbit during the solution iteration process, and the derived along-track 
correction can be expressed as an orbital time-bias for improving subsequent predictions. 

During the prediction process an orbit in the form of geocentric rectangular satellite co- 
ordinates and velocities at 1-minute time intervals and in units of Mm and Mm/day, is 
computed and stored. For each observational epoch, we use 8th order Lagrangean interpola- 
tion to compute the instantaneous position r and velocity v of the satellite. Let r = (x,y,z), 
and v = (v x ,v y , v*), and let the magnitudes of these vectors be r and v respectively. The 
reference frame used is different depending upon the source of the orbit being considered; for 
an orbit generated by numerical integration and using IRVs, the coordinates are given with 
respect to the true equator and mean equinox at Ohrs UT of the epoch of the IRVs. For an 
analytically-derived orbit, the coordinates are given with respect to the mean equator and 
equinox of date. At the observational epoch, we compute the coordinates of the observing 
station in the same reference frame. Let the station coordinates be s = (x s ,y 8 , z 8 ). Then 
the predicted topocentric coordinates of the satellite are 

x T = x - x s , y T — V ~ Vs, z T = z - z 8 , 

and the predicted range is the magnitude of this topocentric vector. We thus form the 
difference of the observed from the predicted range, and in order to solve for corrections to 
the predicted orbit we require partial derivatives of range with respect to the along-track, 
across-track and radial components of the displacement of the predicted orbit. 

2.2 Partial Derivatives. 

(a) Let the along-track displacement or time-bias be AT, 

Then the displaced range R is given by 

R 2 = (x + v x AT - x s ) 2 + (y + v y AT - y s ) 2 + (z + v z AT — z 8 ) 2 

Then 


dAT 


= v x (x + v x AT - x a ) + v y (y + v y AT - y s ) + v z (z + v z AT - z„), 


or approximately 


I 
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dR 

dAT 


= {v x XT + v vVT + v z z T )/R 


(b) Let the across-track displacement be A Q 
The across-track direction is 


r x v 
r • v 


= (/, m, n) , say. 


So, 


l = ( yv z - zv y )frv 

m = ( zv x - xv z )/rv 

n = ( xv y — yv z )/rv 
So the displaced range R is given by 

R 2 = (x + IAQ - x s ) 2 + {y + mAQ - y a f + ( z + nZi< ? _ 2 *)' 


So, 


R- 


dR 


dAQ 

or approximately 


= l(x + IAQ - x s ) + m(y + mAQ - y s ) + + n ^Q z *))’ 


dR 

dAQ 


= (/ xt + m vt + n z t)I R 


(c) Let the radial displacement be Ar 

So the displaced coordinates of the satellite are x + xAr/r , etc. 
So the displaced range R is given by 


n , xAv >,2 i / i ^ ' 

R 2 = (x + x s ) +{y + — 

v r r 


yAr x2 , / , z ^ r v \2 

y a ) + (z + - z t ) 

r r 

Then as before we calculate the approximate partial derivative 


dR 

dAt 


= (xxt + yyr + zzT)/( r R) 
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2.3 Solutions. 

In practice it was found that the across-track correction was always highly correlated with 
the radial correction, and in many cases the solution was indeterminate. Hence both could 
not be solved-for, and so it was decided to suppress the solution for across-track, and solve 
for a radial correction only, which would thus absorb the across-track error. It was also 
found that simple constant along-track and radial corrections to the predicted orbit did not 
m general absorb all the error in the orbits, and that the set of parameters to be determined 
for a particular pass should be selected from along-track and radial displacements and their 
time rates of change and accelerations. We denote these parameters T, T, T, R, R, R. The 
partial derivatives of range with respect to the rates of change and accelerations of these 
parameters were formed from those of the constant terms by mutiplication by t and < 2 , where 
time < is the epoch of each observation relative to the mid-time of the pass. Such a definition 
for the origin of t is optimum in reducing correlations between the unknowns. It was usually 
found necessary to solve for only 4 of the possible six parameters, but occasionally all six 
were required to remove all trends from the observational residuals. 

Initially a scheme was devised to automatically check for very high (>0.999) correlations 
among the 4 unknowns T,T,R, R , and to supress any one or two of them in order to obtain a 
determinate solution, consistent with obtaining a flat track of residuals. However experience 
showed that the values determined for these 4 parameters were always quite small, of the 
order of 0.1 ms/minute for T and 1.0 cm/minute for R, and of similar magnitude for the 
accelerations. We thus imposed a-priori standard errors of these magnitudes upon the 4 
parameters, and allowed the program to solve for all six unknowns. 

In order to iterate the solution, at each stage we replace the predicted coordinates of the 

satellite by the displaced coordinates as determined from the previous solution. So for 
example 


x = x -(- v x AT + IAQ + xAr/r 

and similarly for y and z. It was found that 4 or 5 iterations were usually sufficient, where 
outliers of magnitude greater than 3 x rms were removed at each stage. 

A selection pf some plots of observational residuals is given in Figure 2, where it is clear that 
all trends have been removed; these residuals are readily used to form normal points using 
the Herstmonceux algorithm. - 

As a final check that the track of residuals is indeed statistically flat, we have introduced 
a single-factor analysis of variance test on the residuals. This test checks for significant 
differences between the means of residuals grouped into normal point bins; any differences 
indicate that not all trends have been removed from the residuals, and a warning is given 
that normal points should not be formed until the cause of the problem is traced. It is 

anticipated that this test may prove useful for the detection of system calibration changes 
during a pass. 6 
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3. Conclusion. 

The scheme was found to be determinate for a wide range of satellites and pass durations 
and quality of predictions, and has been adopted at Herstmonceux. An earlier version of the 
software is also in use at several European SLR stations. The deduced values of time bias 
axe used to good effect for improving subsequent predictions, and additionally the software 
has been used to generate time bias values from the quicklook observations from other SLR 
stations, as described by Wood and Gibbs (these proceedings). 
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Figure 2. flange residuals from a pass of ERS-1 (a) and Lageos (b). 
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